M=15;
point = [7, 8];
Sgrid = zeros(M,M);
for i=1:M
    for j=1:M
        Sgrid(i,j) = exp((-(norm([i,j]-point))^2)/M);
    end
end
surf(Sgrid);
[Rgrid, Bgrid] = Bivariate_4point(Sgrid,1/16,3);
n = size(Bgrid,1);
figure(2);
surf(Rgrid);
axis([-1 n -1 n 0 1]);
title('No Boundary Case');
figure(3);
surf(Bgrid);
axis([-1 n -1 n 0 1]);
title('With Boundary Case');